      SUBROUTINE SETUP
      RETURN
      END
C   HS68
      BLOCK DATA
      INCLUDE 'BLOCO.INC'
      INTEGER I,J
      DATA NAME/'HS68'/
      DATA (X(I),I=1,4)/ 4*1.0D0/
      DATA N/4/ , NH/2/ , NG/ 8/
      DATA DEL0/0.05D0/ ,TAU0/1.D0/ ,TAU/.1D0/
      DATA (GUNIT(1,J),J=3,10)/8*1/,
     F     (GUNIT(3,I),I=3,10)/1000,3*1,-1000,3*-1/,
     F     (GUNIT(2,I),I=3,10)/1,2,3,4,
     F     1,2,3,4/
     F     ,((GUNIT(I,J),I=1,3),J=0,2)/-1,0,0,-1,0,0,-1,0,0/
      END
      SUBROUTINE EF(X,FX)
      INCLUDE 'FUCO.INC'
      DOUBLE PRECISION X(NX),FX,A,B,E
      SAVE A,B
      DATA A/2.4D-3/, B/1.D0/
      ICF=ICF+1
      E=DEXP(X(1))-1.D0
      FX=(A-(B*E-X(3))*X(4)/(E+X(4)))/X(1)
      RETURN
      END
      SUBROUTINE EGRADF(X,GRADF)
      INCLUDE 'FUCO.INC'
      DOUBLE PRECISION X(NX),GRADF(NX),A,B,E1,E,Z,D
      SAVE A,B
      DATA A/2.4D-3/, B/1.D0/
      ICGF=ICGF+1
      E1=DEXP(X(1))
      E =E1-1.D0
      Z=B*E-X(3)
      D=E+X(4)
      GRADF(1)=-(A-Z*X(4)/D)/X(1)**2 - X(4)*E1/X(1)*(B*X(4)+X(3))/D**2
      GRADF(2)=0.D0
      GRADF(3)=X(4)/(X(1)*D)
      GRADF(4)=Z/(X(1)*D)*(X(4)/D-1.D0)
      RETURN
      END
      SUBROUTINE EH(I,X,HXI)
      INCLUDE 'FUCO.INC'
      INTEGER I
      DOUBLE PRECISION X(NX),HXI,C,Y,Y1,Y2
      CRES(I)=CRES(I)+1
        C=DSQRT(24.D0)
        IF(I .GT. 1)        GOTO 100
        CALL MDNORD(-X(2),Y)
       HXI=X(3)-2.D0*Y
       RETURN
  100  CONTINUE
       CALL MDNORD(-X(2)+C,Y1)
       CALL MDNORD(-X(2)-C,Y2)
       HXI=X(4)-Y1-Y2
      RETURN
      END
      SUBROUTINE EGRADH(I,X,GRADHI)
      INCLUDE 'FUCO.INC'
      INTEGER I
      DOUBLE PRECISION X(NX),GRADHI(NX),D,C
      CGRES(I)=CGRES(I)+1
       D=DSQRT(6.283185307179586D0)
       IF(I .GT. 1)      GOTO 100
       GRADHI(1)=0.D0
       GRADHI(4)=0.D0
       GRADHI(3)=1.D0
       GRADHI(2)=2.D0*DEXP(-X(2)**2/2.D0)/D
      RETURN
  100 CONTINUE
      C=DSQRT(24.D0)
      GRADHI(1)=0.D0
      GRADHI(3)=0.D0
      GRADHI(4)=1.D0
      GRADHI(2)=(DEXP(-(-X(2)+C)**2/2.D0)+DEXP(-(X(2)+C)**2/2.D0))/D
      RETURN
      END
      SUBROUTINE EG(I,X,GXI)
      INCLUDE 'FUCO.INC'
      INTEGER I
      DOUBLE PRECISION X(NX),UG(4),OG(4),D(4),GXI
      SAVE UG,OG,D
      DATA UG/1.D-1,3*0.D0/, OG/1.D5,1.D2,2*2.D0/, D/1.D3,3*1.D0/
      IF ( GUNIT(1,I+NH) .EQ. -1 ) CRES(I+NH)=CRES(I+NH)+1
      IF(I .GT. 4)       GXI=OG(I-4)-X(I-4)*D(I-4)
      IF(I .LE. 4)       GXI= X(I)*D(I) -UG(I)
      END
      SUBROUTINE EGRADG(I,X,GRADGI)
      INCLUDE 'FUCO.INC'
      INTEGER I,J
      DOUBLE PRECISION X(NX) ,GRADGI(NX),D(4)
      SAVE D
      DATA D/1.D3,3*1.D0/
      DO     100     J=1,4
      GRADGI(J)=0.D0
 100  CONTINUE
      IF(I .LE. 4)   GRADGI(I)=D(I)
      IF(I .GT. 4)   GRADGI(I-4)=-D(I-4)
      RETURN
      END
      SUBROUTINE MDNORD(A,B)
      DOUBLE PRECISION A,B,NORINT
      B=NORINT(A)
      RETURN
      END
C*******************************************************************
      DOUBLE PRECISION FUNCTION NORINT(X)
      IMPLICIT NONE
      DOUBLE PRECISION X
C********* COMPUTES THE GAUSSIAN NORMAL DISTRIBUTION INTEGRAL
C********* PRECISION ABOUT 16 DIGITS
      DOUBLE PRECISION P1(0:8),Q1(0:9),P2(0:5),Q2(0:6)
      DOUBLE PRECISION SQRT2,RWPI,ARG,ARG2,XABS,ERF,ERFC
      SAVE P1,Q1,P2,Q2,SQRT2,RWPI
C********** COEFFICIENTS TAKEN FROM HART COMPUTER APPROXIMATIONS
C********** NO ERFC 5708 AND ERFC 5725
      DATA P1/.37235079815548067D4, .71136632469540499D4,
     F        .67582169641104859D4, .40322670108300497D4,
     F        .16317602687537147D4, .45626145870609263D3,
     F        .86082762211948595D2, .10064858974909542D2,
     F        .56418958676181361/
      DATA Q1/.37235079815548065D4, .11315192081854405D5,
     F        .15802535999402043D5, .13349346561284457D5,
     F        .75424795101934758D4, .29680049014823087D4,
     F        .81762238630454408D3, .15307771075036222D3,
     F        .17839498439139557D2, .1D1/
      DATA P2/.29788656263939929D1, .74097406059647418D1,
     F        .61602098531096305D1, .50190497267842675D1,
     F        .12753666447299660D1, .56418958354775507/
      DATA Q2/.33690752069827528D1, .96089653271927879D1,
     F        .17081440747466004D2, .12048951927855129D2,
     F        .93960340162350542D1, .22605285207673270D1, .1D1/
      DATA SQRT2/1.41421356237390505D0/
      DATA RWPI/.56418958354775629D0/
      XABS=DABS(X)
      IF ( XABS .GT. .5D0 ) THEN
        IF ( XABS .GT. 8.D0 ) THEN
C******* TAKE APPROXIMATION 5725
          IF ( XABS .GT. 100.D0 ) THEN
            ERFC=0.D0
          ELSE
            ARG=XABS/SQRT2
            ERFC=(((((P2(5)*ARG+P2(4))*ARG+P2(3))*ARG+P2(2))*ARG+P2(1)
     F             )*ARG+P2(0))/
     F           ((((((ARG+Q2(5))*ARG+Q2(4))*ARG+Q2(3))*ARG+Q2(2))*ARG
     F                 +Q2(1))*ARG+Q2(0))*DEXP(-ARG**2)
          ENDIF
        ELSE
C******** TAKE APPROXIMATION 5708
          ARG=XABS/SQRT2
          ERFC=((((((((P1(8)*ARG+P1(7))*ARG+P1(6))*ARG+P1(5))*ARG
     F                 +P1(4))*ARG+P1(3))*ARG+P1(2))*ARG+P1(1))*ARG
     F                 +P1(0))/
     F         (((((((((ARG+Q1(8))*ARG+Q1(7))*ARG+Q1(6))*ARG+Q1(5))
     F                 *ARG+Q1(4))*ARG+Q1(3))*ARG+Q1(2))*ARG+Q1(1))
     F                 *ARG+Q1(0))*DEXP(-ARG**2)
        ENDIF
        IF ( X .LT. 0.D0 ) THEN
          NORINT=ERFC*.5D0
          RETURN
        ELSE
          NORINT=(2.D0-ERFC)*.5D0
          RETURN
        ENDIF
      ELSE
C************ TAKE POWER SERIES APPROXIMATION
        ARG=XABS/SQRT2
        ARG2=ARG**2
        ERF=ARG*2.D0*RWPI*((((((((((ARG2/210.D0-1.D0/19.D0)*ARG2
     F     /9.D0+1.D0/17.D0)*ARG2/8.D0-1.D0/15.D0)*ARG2/7.D0+1.D0/13.D0)
     F     *ARG2/6.D0-1.D0/11.D0)*ARG2/5.D0+1.D0/9.D0)*ARG2/4.D0
     F     -1.D0/7.D0)*ARG2/3.D0+1.D0/5.D0)*ARG2/2.D0-1.D0/3.D0)*ARG2
     F     +1.D0)
        IF ( X .GE. 0.D0 ) THEN
          NORINT=(1.D0+ERF)*.5D0
          RETURN
        ELSE
          NORINT=(1.D0-ERF)*.5D0
          RETURN
        ENDIF
      ENDIF
      END


